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SYNOPSIS 


Maleic anhydride is an industrially important chemical. 
It can be produced from oxidation of benzene. Now-a-days its 
production from C^- hydrocarbon has gained considerable 
significance. The most widely used catalyst for this reaction 
is Vanadium- Phosphorous Oxides (VPO) . 

In the present work the kinetics of the reaction in an 
integral tubular reactor, data of which are available in litera- 
ture, has been modelled. Various studies have shown that the 
reaction takes place via a redox mechanism. A model based on 
redox mechanism has been investigated. Using the experimental 
data of Gopichand (1) various kinetic parameters involved in the 
reaction model has been estimated and the plausibility of> the 
proposed model has been discussed. A few other models (1), (4) 
have also been used to estimate the various kinetic parameters . 
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chapter 1 


INTRODUCTION 


Maleic anhyclrjde and its derivatives ate produced 
uorld-wide. This industrially important chemical is used in 
the production of agricultural chemicals, resins, food additives, 
lubricating oil additives and pharmaceuticals. 

Earlier benzene was oxidized to produce maleic 
anhydride. The reaction can be represented by 


C 


6^6 



H - C — C 

\\ \ ^^'2^ + 2C02 

H - C - CC 


( 1 , 1 ) 


Now-a-days benzene is being replaced by -hydro carbon as the 

feed stock. n-Butane feed is generally preferred because of its 

lo^'? cost and further more it is environmentally more desirable 

raw matericil. The catalytic oxidation of butane produces a 

of 

product stream which mainly contains maleic anhydride and carbon 

A 

oxides. The reaction can be represented by 


H 

I 

H- C 
I 

H 


C 

I 

H 


f 

C 

I 

H 


H 

I 

C 

I 

H 


- H t 


0, 


H -C- Cf 

II N 


H -C- C 


0 

0 

/* 

^0 


+ 4H20 


( 1 . 2 ) 


Vanadium Phosphorous Oxide (VPO) system ds generally 
used as the catalyst for the production of maleic anhydride 
from n-butane. In this ca ilyst, many phases may exist, for 
example. Vanadium may show oxidation states li]<e V(III),V (IV), 

V (V) etc» The presence of various phases depend on many 
synthesis variables. The catalytic composition and method of 
preparation were found to be important parameters for the 
reaction . 

In order to understand a catalytic reaction, kinetic 
modelling becomes an important task. In simple terms, modelling 
is to build a mathematical model for the system through which 
it becomes possible to predict, control and optimize the sysrem. 
The modelling of a system is a very involved task. It starts 
with assuming certain plausible models and then estimating the 
parameters appearing in the mov-u i with the available experimental 
data. Next step ±3 to discriminate among the rival models and 
finally check the adequacy of the model thus proposed. 

For the present study, the experimental data were taken 
from literature (1) . These data were taken for a continuous 
reactor. The feed stream consisted of n-butane and air which 
after compression and preheating was passed to the reactor. The 
conversion and yields were large enough, bo that an integral 
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analysivS was possible experiments were conducted in sucVi a way 
that external and intrapellet mass transfer and heat transfer 
resistances were negligible. This is important because in 
order to determine kjnetic parameters, data are required which 
are unmasked by transport resa stances* 

The Vanadium -Phosphorous Oxide (VPO) catalyst used 
by Gopichand (1) had the following composition: 

0 ^ = 35 , 3 % 

O 5 - 9.7 % 

Fe = 0 . 5% 

Silica gel = 74.8% 

The aim of the present work was to model the selective 
oxidation of n-butane to maleic anhydride using Gopichand 's data. 
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CHAPTER 2 


STUDIES OH THE KINETICS AND MECHANISM OF - HYDRO CARBON 
ON VPO CATALYST 


2.1 INTRODUCTION 


The kinetics of the n-butane oxidabion on v.^nadlum 
phosphorous mixed oxides has been studied by Escardino et-al,(2)^ 
and by Wohlfahrt and Hofman (3), among others. The reaction 
scheme proposed was 


n- 


butane 





maleic anhydride 



( 2 . 1 ) 


It was found that at a low hydrocarbon partial pressure, 
the behaviour of n-butane oxidation may be described by three single 
pseudo-first order rate expressions. Oxygen chemisorption was 
found as the rate determining step at large butane pressures. 

Wohlfart and Hofmann studied the reaction using a catalyst ’ 

I 

with a P/V ratio of 1.06, and extended the range of hydrocarbon and ^ 
oxygen concentrations investigated by Escardino et-al (2) and found 


5 


that n~butane is oxidized according to a rate equation of the type 


^ Pb Pq 

L 1- k p 


0,235 


( 2 . 2 ) 


B 


and maleic anhydride is oxidized according to a rate equacion 
of the type 


PmA Po 


0.285 


r = 
MA 


1 -t- k p 


(2.3) 


B 


Where the subscripts B, o and MA de. sign ate butane, 
oxygen and maleic anhydride respectively. 

Another study by centi-et-al (4) gives a different 
reaction scheme. This is given by 


n-Butane 


CO 


maleic anhydride 


CO, CO^ 


(2.4) 


Rate equations of the follov^ing tYP® were proposed by 
centi-et-a] . (4) . 


"" 


^2 = 


Pb Pq 
( 1 + Pb> 


2 -^o 


(2.5) 


6 


e 

-Is 

^3 ^^3 ^MA b 

Pb 

Where is the rate of formation of maleic anhy^^tade 

r^ is the rate of forjnation of Carbon dioxide and is the rate 
of the decomposition of the ma]eic anhydride, 

A redox mechanism has been proposed by Mars and Van 
Krevelen (5). Varma and Saraf (6) aJ so used a redox model to 
fit their data for the oxidation of n~hutcine to maleic anhydride 
on e VPO catalyst. Butadiene was found to be a major intermediate 
is the proposed reaction scheme- 

Butane ^ Butadiene ^ Maleic anhydride 



Carbon oxides (2.6) 

pseudo first order rate expressions were fitted to each 
of the above reaction steps at low hydrocarbon concentrations. 

A study by Cavarii-et-'al (7) has shown that 1-butane 
adsorbed on a V(lV) site oxidizes to butadiene accompanying the 
formation of V(III) site. If the concentration of V(v) sites 
around the adsorbed butadiene is high enough, butadiene is further 
oxidized to maleic anhydride or carbon oxides. On the other hand 
if the concentration of V(V) sites is low, butadiene desorbs from 
the V(III) site. 



7 


In the case of n-butane, when the concentration of n-butane 
is low, intermediates such as butene and butadiene v/ere not detected 
Only at high concentration of n-butane, butene and butadiene were 
detected . 


A redox model for the oxidation of n-butane has been 
profx^sed by Gopichand (1). He proposed that when n-tautane adsorbs 
on V(lV) site, it^onvert^ to olefins and vClIl) forms. If the 
density of the oxidized sites V(v) in the vicinity of the adsorbed 
olefins is high# these olefins are further oxidized to maleic 
anhydride. But if the density of the V(v) sites is low, the olefin.s 
desorb* 


The proposed rate expressions are 


B 


P 02 + Pb 


( 1 + 


'"1 Po 
H Pb 


(2.7) 


2 + ^3 Pb 


. ^"4 PmA, 

k, p ' 

5 


‘MA 


P02 - k 4 


k- p 
1 


^3 Pb 


( 2 . 8 ) 


(1 + 


P, 


MA 


k^ p 


B 


P, 


o. 


b o, 


Pb ^4 Pma 


P 

1 ^o< 


(1 + 


^3 Pb 
’" s Po. 


, ^^4 Pma, 
kT^ 

5 


(2.9) 


^2 Pb 



a 


vJhere r is the rate of reaction of n-butane , 

13 

is the rate of production of oxides of carbon and 
^MA rate of production of maleic anhydride. 

Gopichand (1) obtained the daita using i v/o VPO catalysts 

with and v/j thout iron promoter , using a P/V atomic ratio of 1,20. 

2 / 

The catalyst had a surface area of b 8 m /grn. The data for promoted 

o 

catalyst in the temperature range of 350 C - 400 C were used in 
the present investigation for modelling. The data used are given 
in Appendix 1 . 

2 . 2 SOME GENERAL CONSIDERATIONS 


The analysis is possible either by an Integral method or 
by a differential method. Since the data were available for an 
integral reactor (i.e. high conversion and yields), an integral 
analysis had been used. 

The details of the kinetic model proposed by Gopichand 
as well as other models, tested in the present investigation are 
given in Section 2.3 and 2,4. 

The reactions involved in the selective oxidation of 
n-butane are as follows- 

C 4 + 3.5 O 2 J C 4 H 2 O 3 + 4 H 2 O 

(B) (0) (MA) 


( 2 . 10 ) 


q 


C 4 + 6,5 O. 

(B) (0) 


4 CO 2 + 5 H^O 
(c) 


( 2 . 11 ) 


C 4 + 4.5 0, 

(B) (0) 


4 CO H- 5 HpO 


(c) 


( 2 . 12 ) 


The experimental conditions are such that the ratio of 
the rt 3 te of formation of CO 2 to the rate of formation of CO is 
1,942. 

Hence on the average the reaction for carbon oxide formatior 
can be written as 

C 4 t 5.82 O 2 > 2.64 CO 2 + 1.36 CO t 5 H 2 O (2.13) 

This equation can be obtained from Equations (2.11) and (2,12). 

The partial pressure of butane, maleic anhydride and 
oxides of carbon (p^^/ expressed in terms of 

I 

conversion and yields as follov/s- 

Conversion X of n-butane is given by 

I 

X - _ mol of n-butane converted 

” mol of n-butane fed 1 

^Bo 


Pb = Pbo ^ ^ 


(2.14) 
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of maleic anhyc^nrle 


y 


rnol of maJeic anhydride formed 
mol of n-butane fed 

%A 


so p 


MA 


Pbo y 


(2. IS) 


yield ol oxides of carbon 


z = 


mol of oxides of carbo n i ormed 
mol of butane fed 


P, 


4p 


Bo 


so = (4) (z) p 


Bo 


(2.16) 


Equations (2.10) and (2,13) can be used to find partial 
pressure of oxygen. Let at any point in the reactor, partial 
pressure of maleic anhydride and oxides of carbon be p^^^ and p^. 

Then due to formation of maleic anhydride, the partial 

pressure of oxygen reduces by 3.5 (from Equation (2,10)) and 

due to formation of oxides of carbon, the partial pressure of O 2 
5 82 

reduces by (from Equation (2.13)) so the partial pressure 


11 


of oxygen ccm be \ 7 ritcen as 



P. 


o 


- 3,5 p. 


.82 


MA 


4,0 


P. 


(2.37) 


Where is partial pressure of o'vygen jn feed stream. 

Now on substituting p and p^ (from Equations (2,15) and 

jyiM Cl 

(2,16))in Equation (2.17), we get partial pre.SvSure of oxygen in 
terms of oonversion and yie]ds as 


- p° ~3.5p_ y-5,82 p z (2.18) 

2 ^2 Bo Bo 

Equations (2.14), (2.15), (2,16) and (2.18), assumes 

that there is no volume change during the reaction. This is valid 
since the partial pressure of n-butane in the feed is very low and 
a large amount of inert (nitrogen) is present in the feed. 

We can also v;rite continuity equations for n-butane, maleic 
anhydride and carbon oxides for the isothermal fixed bed reactor. ! 

Consider an element of length 'dl ' in the reactor 

Input of butane = 

Output of butane = (F_ -f dF^) 

y 8 ' 

Disappearance of butane by reaction 

= (-rg) (dW) 

Where (-r_) the mol of butane reacted per sec. per unit mass of the ■ 

I 

catalyst, 1 
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Wow the continuity equation is 

Input = Output h disappearance by reaction -b accumulation 
at steady state, there is no accumulation 
Hence = (Fg i- dF^) 4- i-r^) (dW) 

dFg + (-Ug) (dW) ^ 0 (2.19) 

Now conversion is defined by the equation 

SO dFg = d ( 1 -x)):^ - dx 

using this equation in (2,19) 

-F d X + (~r„) dW = 0 
Bo B 

F_ dx = (-r^) dW 
Bo B 


dx / . 

d“(W/F ' ) ~ 
h)0 


(2 . 20 ) 


Similarly for maleic anhydride formation, continuit', equation is 




= r. 


acw/Fg^) 


( 2 . 21 ) 


and for oxides of carbon, the continuity equation is 


dz 


= r 


cKW/Jbo) 


( 2 . 22 ) 
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In order to solve the continuity equations (?*20), (2,?1) 

and (2,22) we need rate expressions. These rate expressions have 
been obtained from varj ous models proposed for the reaction. The 
solution of these equations give values of conversion and yields 
and these can be used to obtain various parameters involved in the 
react >r>n kinetics. Parameter estiination involves minimization of 
residual sum of squares. 

The initial conditions needed to solve the equations are 


X = 0 
y = 0 
z = 0 


at 



0 


(2,22a) 


2 , 3 MODEL, 1 (REDOX MODEL ) 

The selective oxidation reactions on oxide catalysts are 
assumed to follow redox mechanism. Gopichand (1) has proposed a 
kinetic model based on redox concept involving three types of active 
catalytic sites with different oxidation states, ' 

Vario\is steps in the redox mechanism are as follows- 


n02 

+ 

R - 

-^1 


^2 ^ 

(2.23) 

mB 

+ 

O 2 X - 

^2 


MA + R 

(2.24) 

r B 

+ 

R 

^3 


C t F 

(2.25) 




k , 




SMA 

+ 

R — 

4 


C + F 

(2.26) 




^5 




n02 

-f 

F — 


— > 

R 

(2.27) 
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Where 

B is bubane 

FiA is maleic anhydride 

C is oxides of cax'bon 

R is reduced site 

F .1 h further reduced site and 

^ 2 ^ oxidi'xed site 


n, m, r and s are orders of reaction v/i ch respect to 
oxygen/ butane and maleic anlydride in Ecpaa cions (2.23), (2.24), 

(2.26) and (2,27) respectively. 

Various overall reactions are 


^4 ^10 


C 4 O 3 -I- 4 H 2 O 


(B) 


(MA) 


(2.28) 


C 4 + 5,82 


-> 2.64 CO 2 t 1.36 CO t 5HpO 


(B) 


(2.29) 


C4 H2 O3 + 3O3 


4CO3 t H^O 


(MA) 


(2.30) 


At steady state, the net rates of production of sites 
O^X, R and F are equal to zero. 

From Equations (2.23), (2,24) and (2,28), we can write- 


Oxyqen balance on site X 


’"l PO 2 ®R = *^2 P' 


B 


(2.31) 


Similarly/ oxygen balance on site R(u.sjng "Equations 2.23 t.o 2.30) 
Pqj ®R ^3 P|3 ®R ^4 PhA ®R 


= K '■ P 02 < ’ - ®R - ®02X> 


Using (2.31) in the above ecxuation 


? ’^3 Pfl ^R"* ^ ^4 PmA ®R '^5 Pog ^^"®R~ 


(2.32) 


put the value o£ 0 from (2.31) in Equation (2.32) 

o ^ ^ 

T n 

p 0 

? ^3 ®R ^4 %A ®R ^3 ^Ory ^5 JZ 

2 2 c<lC2 Pb 


®r(P’"3 Pb \ +^’"4 PmA ’^5 pS 


’'5 Pq^ Po.,\ 

+ V pH 

2 ,< / 5 02 


®R = 


k. p" 

5 o, 


(■p kj pP +>? k^ H- kg + kg p"^ 


k. p’^ 

£ 

, , m 

^kjPB 
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Poj Pb ^ ^^4 Pma 

1 n ^ , n^ 

^’^2Pb ’S Po. ^5 Po. 


(’.33) 


In the nbove equations, C< / P and 9 ere stoichiometric 
coo rf 1 cients . From Equation (?,28)o(^is equal to 3,5, from 
Equation (2.29) p is equal to 5.82 and from 3iC[uaT:ion (2.3 0) ^ 
is equal to 3.0, 


Now, from Equations (2.24) and (2.25), the rate of dcpietlon 
of n-buLane can be written as 


m r 

^ ^2 ^3 ®R 


Using Equation (2,31) for 0 

, n 2 

~ ^2 T~ m ®R 


■2 


^ ^2^B 


■3 '"R 


n 


Pq. 




+ 


Using the vaiue of 9^ fr'Cit Eqn, (2.33) in the above 

r\ 

expression, we get 


Pb = ^ - 


*■ ^3 Pb^ 


P. 


(it 


■P^3 P 


^ kop' 


m 
2^B 


IT 


^"5 Pq. 


^ ’"4 Pma 

^5 Pq^ “ 


]7 


or 


-r 


D 


k- 

( ] + — 
^2 

Po 

2 + 
m 

Pb 

VJhere k| 

^1 

< 


= 13 k. 

*^4 

= k 


n 


r 

^3 Pb 


^^3 Pg 

5 ^ o. 


^^4 ^riA 
5 


(2.34) 


From Equcitions (2,24) and (2,26), the rate of formation 
of maleic anhydride can be written as 


^MA ^2 ®025C ” ^^4 ^14 A ®R 


Using Equation (2 .31) for 0 


02 X 


- n 

= k, Po 60 - k. p.® Sn 

2 / , m R 4 ^MA R 

MA ^ ^2 Pg 


= ( 


Po 


2 . S V ^ 


using Equation (2.33) for 0 


R 


J8 


= ( 


1 o. 


'^4 




Po„ 

( 1 4 ± + + 

’^2Pb ’^sPoj 


^^^4 Pma 

— -~ir> 

’^5 Po. 


bo, 


HA 


P 


1 P 


^4 P 


MA 


(1 4 


■» 

^2Pb 


>^^3 Pb 

5^0. 


^"^4 Pma 

Po. 


(2.35) 


From Equations (2.25) and (2,26), the rate of formation 
of oxides of carbon can be written as 


^3 Pb ®r '*' ^4 Pma ®r 


Using Equation (2,33) for 6 


R 


(ka pg I- 


Po2 ^ >^3 Pb 9 »'4Pma i 

( 3. + 1- — — + — ^ ) 

■I n 

kc- P 
5 ^Oo 


y , m 1 n 

^"2Pb ^bPq. 


^3 Pb ^4 Pma 


(i + 


k: ' 

Pq. 


k^ p| 


m 


’^2 Pb 


n 


>"6 Po. 


>^4 Pma 


n 


Po, 


(2.36) 
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The rate expressions (2,34), (2,35) and (2.36) can be 

expx-essed in terms of conversion of n~butane X, yield of maleic 
anhydride y and yield of carbon oxides z, using Equations (2.14), 
(2.15) and (2,18). 

The resultinq rate expressions in terms of conversion and 
yields are used ro solve continuity Equations (2,20), (2.21) and 

( 2 . 22 ) . 

Gopichand (1) assumed the reaction orders n, rn, r 
and s involved in Equations (2.23), (2.24), (2,25) and (2.26) 

equal to one. Using thj s assumption, the rate expressions (2,34), 
(2,35) and (2.36) become same as proposed by Gopichand (l) and 
given by Equations (2,7), (2,8) and (2,9) 

In the present investigation, parameter estimation has 
been done using Gopichand 's model as vjell as above proposed 
general model which involves various reaction orders. 

Since the oxygen present in the feed stream is in abundance 
the reaction order with respect to oxygen 'n' has been taken as zero 

Also various investigations have shown that the order of 
reaction with respect to butane *m' (Equation 2.24) is one. The 
order ’r' and 's' in Equations (2,25) and (2,26) varies from 0.5 
to 0.25, 

Various combinations of 'r' and 's* have been used in the 
parameter estimation. Details of the results are included in 
Chapter 4. 
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2.4 MODEL 2 (CENTT'S MODEL) 


Cent! et-al (4) proposed a model for the reaction for 
conversion of the order of 10%. The reaction mechanrsm is 

n-Butane ^ Maleic Anhydride 

CO 2 

n- Butane is oxidized along two parallel routes* The first 
route, to maleic anhydride depends on the hydrocarbon concentration 
according to Langmuir Hinshelwood model, whereas the second route, 
to carbon dioxide does not depend on the hydrocarbon concentration 
Increasing the n-butane concentration increases the selectivity in 
maleic Anhydride formation. Both reactions have slight dependence 
on the oxygen concentration, except at very low oxygen values, also 
the selectivity in maleic anhydride is constant in the range of 
experimented O 2 concentrations . Their proposed equations are 



r 


1 


- r 


ma 


^1 


(2.37) 


2 CO2 1^2 


(2.38) 




Bu t th'ese Gqiiations \/ere found to be inadequate when 
kinetic analysis was extended to integral conditions in the 
0-100% conversion range , These equations were found to be 
adequate upto 60% coversion. For higher conversions, the 
experimental yields of maleHc anhydride were lower than the 
calculated ones and for conversion of about 100%, the experi- 
mental yield of maleic anhydride is very lo^/. Therefore a 
success! VG reaction of combustion of cte maleic anhydride must 
occur . Parallel to the decrease in yield of maleic anhydride, 
a strong increase in the yields of both carbon mono and dioxides 
was found. In order to take these facts into account, the 
complete modified reaction model can be expressed as- 


1 3 

n- Butane ^ Maleic Anhydride, >00/ CO^ 


CO, 



rd 

The rate expression for the 3 reaction is 


^3 = = ^"3 PmA 




(Pb) 


(2.39) 


The rate equation of n-butane oxidation is 
r. 


CO, 


•r = r 
B mA 


+ 


4 


(2.40) 
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CHAPTER 3 

PARAMETER EbTJHA^'.TlOM 

3 .1 INTRODUCTION 

The experimental data of Gopichand have been used for 
modelling. The data at three different temperatux-es and various 
partial pressures of butane in the feedstream are given in 
Appendix 1 . The data consisted of outlet yields and conversion 
at various W/F , The contiuni ty Equations (2,20), (2.21) and 
(2,22) are also in terms of conversion and yields. For a given 
set of parameters, the solution of these continuity equations 
give values of conversion and yields at various W/f^ These are 
calculated values based on a particular model. These values on 
subtraction from the experimental values give the deviations of 
calculated values from the experimental values. The purpose is 
to find a set of parameters whi ch minimizes the objective function 

F = E. (x, ^ X, - ^i 

i=l L meas . cal meas cal 

+ (z. - z. )^ 

■*■01635 cal 


(3 ,1) 


21 


"th. * 

Where indicates i^ experiment 

meas. - for exi^erimenta]. value 

calc. - for calculated value 

M - total number of experiments 

This function F is to be minimized v/n th respect to various 

parameters . 

Several methods are available to minimize the function. 

Since these kinetic parf-^eters can only have values greater than 
zero, the constrained minimization technique can also be employed. 
However, for the present study, unconstrained minimization method 
has been used. 

The unconstrained methods can be classified into two broad 
categories as direct search methods and descent methods (8) . 

Generally the descent methods are more efficient compared 
to the direct search methods, but it requires, in addition to 
function evaluations, the evaluation of first and possibly higher 
order derivatives of the objective function F. 

On the other hand, the direct search methods require only 
objective function evaluation and do not use the partial derivatives 


of the function 


P4 


Since in the present case, the exact functional relationshit 
for conversion and yields in terms of kinetic parameters is not 
known, the derivative of the ol'*j active function ^/lih respect to 
parameters ccn. not be calculated. So descent method, although 
desirable, car not be used. 

For thj s reason, direct search method has been used for 
minimization of 'P', 

Among the various direct search methods available, it 
is the simplex method v/hich is most accurate and efficient. So 
a simplex method has been used v^hich is described in the following 

I 

section . 

3.2 SIMPLEX METHOD 

Simplex is defined as a geometric figure formed by a set 
of (nH-1) points in the n-da mensional space. For example, in a 
two dimensional space, the simplex is a triangle. 

The basic idea in the simplex method is to compare the 
values of the objective function at the (n+1) vertices of a 
simplex and move this simplex gradually towards the optimum point 
during the iterative process. The movement of the simplex is 
achieved by using three operations known as reflection, contraction 
and expansion. The details of the method are given in Appendix 2. 
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3.3 SOLUTION OF DIFFERENTIAL EQUATIONS 


The minimization toy implex method require values of the 
otojective function. Since \je have only differential equation for 
the conversion and yields, the value of the objective function can 
toe determined only after solving the equations . 

The set of differential equations can be solved with the 

"f* It. 

given initial conditions by the 4 -order Runga Kutta method. 
Consider a set of 'N' differential equations 


ar = fj ( t, y ) 
dy^ 


(3.2) 


II 


dt 


fj (t, y) 


dt" = 

with initial condition y = 0 at t = 0 
where y = (y^ , ^ 2 * * 

The 4 -order Runga Kutta method gives the values of 
various functions y^, 72 (to+1) point when values of these 
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th 

functions at n point is knovm. The solution is gaven by (3*3) 


yj,n+i = y 


J, n 


+ § (k 
6 


J, 1 ^ 2 ^ ^ J# 3 4^ 


(3.3) 


J 1/ a««>aI''T 


'''('here 


th 


Yj : Value of function at (n+l) point 


th 

Yt ; Value of function y_ at n point 

J / n ij 


k 


J,1 


= • U.n • y2.n ' ••••yN.n’ 


J = 1, 


,N 


'J, 2 


h 


h 


2 ' ^l^n 2 ^1] ' ^2,n 2 ^21 ‘ 


' %,n 2 ^N.n^ 


■J,3 


+ I '• yi,n ' I ’"12 ' ^2,0 


, h T 

+ 2 ^2? ’ 


‘ ' ^N,n 2 ^N, 2^ 


^,4 


= fj(t^ h h ; t h ; y^^^ + h 1^23 ; 


• %,n ^ 


J = 1, 2 , 


.N 


h is the step size 


The fiow chart describing the process of parameter 


estimation is shown in Pig. 3,1. 
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InjticiJ guess Tor the meters 



Fig, 3,1. Flow Chart for Parameter Estimation 
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3 . 4 SOME 0THE:R ASPECTS 


From physical considerations, it is clear that the 
conversjon Oj' Butane in the reactor can not exceed unity. This 
feictor is important and has been taben into account while 
determinang the parameters. 

After each iteration in che minimization process, a new 
set of parameters js generated. The conversion and yield corres- 
ponding to this set of the parameters is obtained by solving the 
continuity equation using Runga Kutta method. If conversion is 
found to be greater than one, the minimization is stopped and the 
value of parameters at that iteration is taken as the final values 
even if it does not satisfy the convergence criterion. 

If the conversion is < 1 and the convergence criterion 
is not satisfied then the iterative process is repeated. 

The computer programme for parameter estimation hris been 


given in Appendix 3 
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CHAPTER 4 


RESULTS AND DISCUSSION 


The aim of the present study was to model the selective 
oxidation of n~butane to maleic anhydride using Gopichand's 
experimental data (Appendix 1) . A general model based on redox 
mechanism has been proposed }Dy Gopichand (1) . This has been 
slightly modified taking into account the orders of the elementary 
reactions in the proposed model. This model is given by Equation 
(2,20) to (2.22); (2,34) to (2,36) and (2,22a). The same equation 
apply to Gopichand's model except that in Gopichand's model, the 
orders of reactions m, n, r, s are all equal to unity. In a third 
model viz., Centi's Model (4) , the rate expressions for the 
productions maleic anhydride , CO 2 and rate of consumption of 
n-butane are given by Equations (2.38) to (2,40) in place of 
(2.34) to (2,36). The results obtained vrith different models 
are discussed in the following sections. Parameters involved in 
this model as well as in the other model proposed by Gopichand (1) 
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arifi Cent.l,.et, al (4) have been esi:irnateci using sequential 
ntinple:: methort. 


^ ST-\R3^IMG SIMPLICES FOR PARAMETER b.S ^^IMATTQN 

For the general moclel as v;ell as Gopichand's model, the 
starting simplices at vorions temperatures (350°C, 375°C and 
400°C) are given in Tables ? and 3 respectively* For Centl's 
model, the starting simplex is given in Table 4. 

d . 2 GOPTCHAND ' S MODEL 


The va].ues oP the intrinsic parameter at different 
temperatures using Gopichand*s model are given in Tabie 5* The 
values of all the parameter were found to be positive* The rate 
constants are plotted in the form of Arrhenius plots, which are 
given in Fig* 4.1* The activation energies are given in Table 5, 
The activation energies for 1^2/ and are rather low. This may 
be due to the nature of data (integral) used, as well as the 
uncertainities associated with nonlinear esbimation. On the 
average each parameter estimation at a given temperature for a 
given model took approximately 30-40 min* on PC - XT. 
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4.3 GENERAL MODEL INVOLVING ORDERS 


VJhtle estimal'ing parciineters orders 'm' and 'n' have 

always been taken as one and zero respectively. Mso, four 
combinations of 'r’ and 's' have been used. The value of 'r' 
and 's' were taken as 0.5 or 0.25. 

The value of kinetic parameters for various combinations 
are given in Table 6 to 9. 

From Table 6 and 7 it can be observed that by changing 

's' from 0.5 to 0.25 and keeping 'r' equal to 0.5# the value of 

parameters did not change. Also from Tables 8 and 9# the value 
of parameters remain unchancjed by changing 's' from 0,5 to 0,25 

while keeping 'r' as constant and equal to 0.25. 

I 

These observations lead to the conclusion that the 
kinetic parameters were not much influenced by the order 's'. ^ 

i 

f 

Tlie rate constants have been plotted in the form of I 

Arrhenius plots (In k vs ~) and these plots are shown in Figures 
4.2 and 4.3. The activation energies E 2 and E^ are relatively 
small . 

i 

* In general# the value of rate constant should increase 

li 

r. 

with temperature. However in table 5-9# some spurious values of ; 

i 

( 

rate constants were obtained# which are physically not possibly. 

; , ' , I 

Be titer data as' well as better estimation techniques are desirable* | 

' '■ ■ { 

c 
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4*4 CENTI ' S MODEL 

The values of Iclnetlc parameters at various 

1 

temperatures are given in Table 10. A plot of In k vs — is 
given in Fig. 4.4. The activation energy corxesponding to 
various parameters is of the order of 10 kcal/mol which is 
very small. So it can be concluded that Gopichand's data do 
not fit well with the Centi's model. 


4.5 CONCUSSIONS 


Prom the above discussion it can be concluded that 
Gopichand's data best fit in the redox model. The centi's model 
based on Langmuir-ad sorption concept is found to be not much 
useful . 

The activation energies as well as R.G.S, values are 
comparable in the Gopichand's model and general redox model. 

So from these data it is not possible to discriminate among 


these two models 
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Table 1 


O 

Initial Simplex Corresponding to Temperature 350 C for Redo:>c Model 



1^2 

’^3 

’^4 


1.31 

0.07 

125.34 

0.43 

1J6.43 

5.92 

0.03 

134.54 

0.39 

115.60 

2.54 

0.63 

295.54 

0.75 

135.65 

4.64 

0.05 

165.84 

0.32 

123 .55 

1 .04 

0.04 

145.84 

0.63 

]64.93 

o 

o 

» 

0.05 

234.54 

0.58 

129.04 
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Table 2 


Initial Simplex Corresponding to Temperature 375° C for 
Redox Model 


k, 

to 

^3 

^4 

^5 

0.1357 

0.Q439 

5 .2439 

0.3268 

21.1234 

0.3542 

0.0456 

7.4328 

0.3216 

29.4231 

0.3543 

0.2356 

9.8065 

0,3218 

34.5264 

0,2467 

0.5086 

7.4321 

0 .5436 

28.4356 

0.4684 

0.0454 

5.4728 

0.7216 

30.4231 

0,3487 

0.9643 

7.5670 

0.6742 

36.1247 




Table 3 


Initial Simplex Corresponding to Temperature 40C C for 
Redox Model 


’a 


'*"3 

^4 


0.03 


0.56 

0.54 

8.97 

0.43 

0.03 

0.34 

0.87 

3.65 

0.32 

0.02 

0.54 

0.21 

4.65 

0.05 

0.04 

0.54 

0.32 

1.54 

0.78 

0.34 

0.09 

0.43 

10.02 

0.45 

0.94 

0.03 

0.54 

4.75 
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Table 6 

Kinetic Parameters for General Redox node!) 


(orders r = 

0.5, s = 

0,5, m = 1.0, 

n = 0) 


Tempera ture 

Parameters 

350°C 

375 C 

400''c 

Activa tion 
Energy 
kcal/mol 

ki 

0.3480 

0.2057 

3 .4320 

34 .38 


0.2780 

0.7570 

0.7657 

17,46 

^3 

0.3120 

7.446 

324.50 

113 .03 


0.5400 

0.5046 

0.7008 

3 .90 

^5 

5.7860 

30.1900 

109.90 

48.46 

R.S.S. 

2.0267 

1,1828 

0.6539 



2 _ R,S.S . 
^ “ n - k 


0.1351 


0.0788 


0.0436 



39 


Table 7 


Kinetic Parameters for General Redox Model 


(orders r 

=0.5, s = 

0 .25 , m = 1.0, 

n = 0) 

Temperature 

Parameters 

350°C 

375°G 

400° C 

ki 

0.3480 

0.2057 

3.432 


0.2780 

0.7570 

0,7657 

^3 

0.3120 

7 .446 

324.5 


0,5400 

0.5046 

0.7008 


5.7860 

30.19 

109.9 

R.S.S. 

2,0267 

1.4770 

0.6539 


2 R«S.S . 
^ n -k 


0.1351 


0.09S5 


0.0436 
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Table 8 

Kinetic Parameters for General Redox Model 
(orders r = 0.25, s = 0.5, m = 1.0, n = 0) 


Tempe rature 
Parame ters 




R ,S .S . 


2.0267 


1.3996 


0.9952 



0.1351 


0.0933 


0.0663 
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Table 9 


Kinetic Parameters for General Redox Plodel 
(orders r = 0.25, s = 0.25, m = 1.0, n - 0) 


Tempern bnre 

Parameters 


375°C 

40o'‘c 

Activation 

Energy 

kcal/rnol 

kj 

0.3^60 

0,3145 

2.201 

29,42 

^2 

0.2760 

0.3509 

0.6778 

14.35 


0.3520 

7 .6420 

193.00 

102.95 

^4 

0.4080 

0.5 166 

0.5712 

5.62 

^5 

5 .9860 

31.7900 

134.20 

51,05 

R.S.S. 

1 .9099 

1.5137 

0.9952 


.2 _ R.S.S. 
n-k 

0.1273 

1.0091 

0.0663 



rENTr^'^L LlBRAi 

1 1 T l^ANPUn 




Acc. No. 
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TabJe 10 


Klnebic Parameters for Centj ' s Model 


Tempera ture 

Parameters 

350°C 

375*^0 

4no”c 


0.7466 

0.5924 

0.9470 


0.7467 

0.8833 

0.5958 

’^2 

1 .3240 

3.1250 

2.7170 

^3 

1 .0460 

0.9169 

1.7281 

a 

0.4270 

0.3167 

0.5138 

b 

0.0735 

0.0073 

0.1968 

c 

4.1820 

4.7060 

4.2190 

d 

1.2200 

1.0500 

1.3800 

R.S.S. 

2.1526 

3.8349 

6.8646 

_ R.S.S. 
n -k 

0.1794 

0.3196 

0,5721 


Activation Energy 


E ^ 3 « 25 

E^ = 12.66 


7.51 


kcal/mol 

kcal/mol 

kcal/mol 
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APPEITDi;; 1 


Table 11 


Experimental Data at Temperature 350°C 


w 

F 

Bo 

PBO* 
a tm , 

Conversion 

of 

B 

Yield 

of 

MA 

Yield 

of 

C 

3.0 

0.006 

0.4317 

0.1205 

0.3166 

4.5 

0.006 

0.5317 

0.1558 

0.3763 

6.0 

0.006 

0.6267 

0.1933 

0.4329 

7.5 

0.006 

0,6900 

0.2183 

0.4625 

9.0 

0.006 

0.7350 

0.2450 

0,4917 

3.0 

0.008 

0.4800 

0.1095 

0.3706 

4.5 

0.008 

0.6113 

0.1500 

0 .4613 

6.0 

0.008 

0.6725 

0.1725 

0 ,5013 

7.5 

0.008 

0,7388 

0.1987 

0.5396 

9.0 

0,008 

0.7863 

0.2200 

0.5672 

3,0 

0.010 

0.5330 

0.1040 

0.4286 

4.5 

0.010 

0.6450 

0-1350 

0.5110 

6.0 

0.010 

0.7190 

0.1580 

0,5605 

7.5 

O.OlO 

0,7710 

0.1780 

0 .5940 

9.0 

0.010 

0.8250 

0.2000 

0 ,6250 

3 .0 

0.012 

0.6292 

0.1100 

0.5185 

4.5 

0.012 

0.6875 

0.1258 

0,5619 

6.0 

0,012 

0.7682 

0.1500 

0 ,6146 

7.5 

0‘.012 

0.8200 

0,1683 

0,6500 

9.0 

0.012 

0.8392 

0.1767 

0,6625 


* PBO is partial pressure of butane in feed stream* 
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Table 12 


Experimental Data at Temperature 375 "c 


w 

^Bo 

PBO* 

atm. 

Conversion 

of 

B. 

Yield 

of 

MA 

Yield 

of 

c 

3.0 

0.006 

0 .2250 

0.332S 

0.0928 

4.5 

0 .006 

0.3017 

0.1800 

0.1212 

6.0 

0 .006 

0 .3883 

0.2367 

0.1513 

7.5 

0.006 

0.4450 

0.2733 

0.1704 

9.0 

0 . 006 

0.5083 

0.3183 

0.1906 

3.0 

0.008 

0.2413 

0.1250 

0.1159 

4.5 

0.008 

0.3488 

0,1850 

0,1634 

6.0 

0 ,008 

0.4088 

0.2200 

0.1884 

7.5 

0.008 

0.4863 

0.2675 

0.2186 

9.0 

0.008 

0.5500 

0.3075 

0.2422 

3.0 

0.010 

0.2650 

0.1230 

0.1420 

4.5 

0.010 

0.3610 

0.1710 

0.1895 

6.0 

0.010 

0 .4400 

0.2140 

0.2265 

7.5 

0.010 

0.5070 

0.2500 

0.2560 

9.0 

0.010 

0.5880 

0.2980 

0,2898 

3.0 

0.012 

0.3300 

0.1400 

0,1900 

4.5 

0.012 

0.3883 

0.1675 

0.2210 

6.0 

0.012 

0.4842 

0.2142 

0.2697 

7.5 

0.012 

0.5608 

0.2542 

0.3065 

9.0 

0.012 

0.5933 

0.2717 

0.3196 

* PBO 

is partial 

pressure of 

butane in 

feed stream. 


Table 13 


Experimental Data at Temperature 400 °G 


w 

H' 

Bo 

PBO * 
atrn. 

Conversion 

of 

B 

Yield 

of 

MA 

Yield 

of 

C 

3,0 

0 . 006 

0.1433 

0.1223 

0.0215 

4.5 

0.006 

0.1983 

0,1667 

0,0316 

6.0 

0.006 

0.2667 

0.2217 

0,0451 

7.5 

0,006 

0.3133 

0,2567 

0,0550 

9.0 

0.006 

0.3700 

0.3017 

0.0683 

3 .0 

0.000 

0.1433 

0,1139 

0.0268 

4.5 

0.008 

0.2150 

0.1712 

0.0441 

6,0 

0.008 

0.2600 

0.2050 

0.0559 

7,5 

0.008 

0.3238 

0.2500 

0.0728 

9.0 

0.008 

0.3813 

0.2900 

0.0900 

3 .0 

0.010 

0.1440 

0,1110 

0.0329 

4.5 

0.010 

0.2070 

0.1570 

0.0504 

6.0 

0.010 

0.2650 

0.1970 

0.0603 

7.5 

0.010 

0.3190 

0.2330 

0-0055 

9.0 

0.010 

0.3920 

0.2800 

0.1115 

3.0 

0.012 

0.1733 

0.1267 

0.0470 

4.5 

0.012 

0,2117 

0.1517 

0.0593 

6.0 

0.012 

0.2817 

0.1975 

0.0839 

7.5 

0.012 

0.3458 

0.2367 

0.1087 

9.0 

0.012 

0.3750 

0,2542 

0.1207 

* PBO 

is partial 

pressure of 

butane in 

feed stream 
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APPENDIX 2 


SIMPLEX METHOD 

In an n dimensional space, a geometrac fi-jnre formed 
by a set of (nH-1) points is called a simplex. In the simplex 
method the values of the objective function at the (n+1) vertices 
is compared and the simplex is moved to an optimum point during 
the iterative process. The movement of simplex involve three 
operations reflection, contraction and expansion, the details 
of which are given by Rao (8) and are as follows: 

(i) REFLECTION 

If is the vertex corresponding to the highest value 
of the objecLive function among the vertices of a simplex, one 
can expect the point obtained by reflecting the point in 
the opposite face to h~ve the smallest value. If this is the 
case, then one can construct a new simplex by rejecting the point 
Xj^ from the simplex and including the new point X^. 

This procedure is shown in the Fig. A(la). The points 
X^,/ form the original simplex, and the point X^, X 2 and X^ 

form the new one. Similarly in Fig.^Adb), the original simplex 
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is given by points X^, end and the new one by X^, 

and X^ , Again one Cim construct a new simplex f rora the 
present one by rejecting the vertex corresponding to the highest 
functional value. 

Since the direction of movement of the simplex is always 
awoy from the worst result, we wil] be moving in a favourable 
direction. If the objective function does not h^''ve step valleys, 
repetitive application of the reflection process leads to a 
xig-zag path in the general direction of the minimuin r:.s shown 
in Fig , A (2) , 

Mathematically, the reflected point X^ is given by 

= ( 1 -1-^) X^ (1) 

'Afhere X|^ is the vertex corresponding to the maximum function 
value i,e, 

f (Xj^) = max f (X^) ? 
i=l to (n+1) 

X is the centroid of all the points X. except i - n and 
o 


X 


n+1 

n 



i=l 

i^h 


given by 


1 

n 
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and <*4, 0 is the reflection coefficient defined as 

distance between X and X 
^ = r o 

distance between Xi and X 

n o 

will lie on the line joining and on the 

far side of X^ from X^^ with |x^ - X^| = c< jx|^ x^| , If f(X^) 

lies between f f(Xj^) where X^ corresponds to the 

minimum function value i,e, 

f (X^ ) = min f (Xj^) 

i~l to (ntl) 

X|^ is replaced by and a new simplex is started. 

If we use only the reflection process for finding the 
minimum# we may encounter certain difficulties in some cases. 

For example if one of the simplexes (triangle in two dimensions) 
straddles a valley as shown in Fig.A(3) and if the reflected 
point X^ happen to have an objective function value equal to 
the of the point X|^# we will enter into a closed cycle of 
operations. Thus if X^ is the worst point in the simples defined 
by the vertices ^ 3 ' reflection process gives a new 

simplex with vertices X^, and Again has the highest 
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fmictionEil vaJue out of the vertices and we obtain 

old simplex itself by using the reflection process. Th\is the 
optimisation process is stranded over the valJey and there is no 
way of moving towards the or'timum point. This trouble can be 

overcome by malcjng a r- thst no return can be made to points 

\ 

which h<ive just been 



Whenever such sj ^xion is encountered, we reject the 
vertex corresponding to the second worst valiie instead of the 
vertex corresponding to the v;orst functio 2 :i value. This method, 
in general, leads the process to continue towards the region of 
the desired rnjn: 2 mum.» However the final simplex may again 
straddle the minimum, or it may lie within a distance of the 
order oE its own size from the minimum. In such cases, it may 
not be possible to obtain a new simplex with vertices closer 
to the mimimum compared to those of the previous simplex, and 
the pattern may le td to a cycle process a sho\Ai in Fjg. A (4). 

In the above, sucessive simplexes formed from the simplex 
123, 234, 245, 456, 467, 478, 348, ... which can be seen to be 
forming a cycle process. Whenever this type of cycling is 
observed, one can talce the vertex that is occuring in every 
simplex (point 4 in fig.) as rhe best approximation to the 
optimum point. If more accuracy is desired, the simplex has to 
be contracted or reduced in size. 


I 
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(ii) LXPANSION 


If a roflec-tion process gives a iJOint for \'Jhich 
^ (V <C f(X 3 _) i.e. if the reflection prod\ices a new minimurriy 
one can generally expect to decrease the function value 
further by moving along the liirectjon and X^. Hence we 
expand to by the relation, 

X^=^X+(l-'i))X (2) 

e r o 

r, is called expansion coefficient defined as 

^ distance between and X_ 

9 ® 9 - > 1 

distance between X and 

r o 

if f (X^) < f(X^), we replace the point X^ by and restart 
the process of reflection. On the other hand if f (X^) > f (X^^) 
it means that the expansion process is not successful and hence 
we replace bhe point by X^ and start the reflection process 
again . 

(iii) CONTRACTION 

If the reflection process gives a point for which 

f(X ) > f(X,) for all i excepting i = h and f (X^) < ^ 

then we replace the point by X^, Thus the new X|^ will be 

X , In this case the simplex* is contracted as follows; 
r 


P \ + (, 1 -?- ) 


( 3 ) 
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\^here P Is called contraction coefficient 04 1 and is 

defined as 

. distance bet;/een X and X 
p c o 

distance between X, and X 

n o 

if ^ ^ f / use Equation (3) \;ithout changing the previous 

point X^ , If the contraction process produces □ point for 
which f ^ rriin [JflXj^), f (X^) 3 , we replace the point 

in , X^, .... by and proceed with the reflection 

process agesin. On the other hand if f( X^) ^ min 
f (X ) J f the contraction process will be a failure and in this 
case, we replace all X^ by (X^ + X^)/2 are restart the 
reflection process. 

The method is assumed to have converged whenever the 
standard deviation of the function at the (n-tl) vertices of 
the current simplex is smaller than srxne prescribed small 
quantity i.e. 


n-f 1 
i=l 


f(X^) f(X^) 


^2 


Q 




Fi g.. A ( 2) Progress of reflection process. 





APPRNDIX 


- ^ 

C PARAMETKR ESTIKATION 


C THIS PROGRANnE SOLVE THE CONTINUITY EQUATIONS USING TOURTH ORDER 

C RUMGA KUTTA METHOD AND USES SIMPLEX METHOD TO MINIMIZES THE LEAST 

C SQUARE FUNCTION TO ESTIMATE THE PARAMETERS. 

C 

REAL Xf20I»PB0f20'),Yn('20,3I 

REAL P(6,5),Y(6),yif3),Kf5),YCf20,3).XCC20I,KCr5),PARf6I 
INTEGER M,NDIM,L 
REAL H. MAX. SUMP 
C 

C ’PBO’ IS PARTIAL PRESSURE OF BUTANE IN FEED STREAM. 

C ’YM’ IS A MATRIX OF CONVERvSION AND YIELDS. 

C X=W/FO UHERB 'U» IS HEIGHT OF CATALYST FO ' IS FLOW RATE OF FEED STRI , 

C 'P' IS A MATRIX OF INITIAL PARAMETERS. 

C ’Y’ IS A VECTOR CONTAINING VALUE OF FUNCTION CORRESPONDING TO INITIAfi 

C PARAMETERS. 

C ’YC’ IS A VECTOR CONTAINING INITIAL CONDITIONS TO SOLVE THE CONTINUn 

C EQUATIONS. 

r 

C 

OPENf UNIT=r> , FILE*-^ ' FOR.'') . DAT ' ) 

OPENCUNIT^6 , FILE« ' FOR6 . DAT » ) 

OPEN(’UNIT=^21 ,FILE'’FOR21 .DAT’ 'I 

M-20 

NDIM-^5 

MP^CNDIM+tl 

NP--NDIM 

FTOL^O.Ol 

DO 20 1=^1, M 

READ(5,*) XCI),PBO(J).(YM(I,J).J-l,3) 

20 CONTINUE 

DO 30 I=-l.(MDIM+l 1 
READ(6, A)fP(I , J), J«1,NDIMI 
30 CONTINUE 

DO 50 1=1,6 
DO 40 J=l,5 
KCJ)=P(I,J) 

40 CONTINUE 

CALL RUNGA(K,YC,XC‘) 

DO 10 L=l .20 
10 CONTINUE 

Yf ii=funk(:yc,yii) 

50 CONTINUE 

C URITEC21,996) 

CALL AMOEBAf P , Y ,MP , NP ,NDIM , FTOL , ITER , YM) 

C URITEr2t,996) 

C DO 60 1=1,6 

C URITEC21,*)(P(I , J) 5) 

C WRITRf 21,998) YCH 


1 
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C 60 


100 


?.00 


400 

300 


500 


600 


C 

996 

997 

998 
1000 
JlOO 
1200 


C 

C 


30 

20 


C 

C 


CONTINUE 

MAX^Yri) 

DO 100 1-2,6 

IFrYfn.GT.MAX') MAX = YCI) 

CONTINUE 

DO 200 J-1,6 

IFCYCJ) .Ep.MAX) 1.^.7 

CONTINUE 

DO 300 1-1,5 

DO 400 J=<1,6 
IF(J.EQ.L)GO TO 400 
SUMP-SUMP-JPCJ, I) 

CONTINUE 

PARC I)= (SUMP/ 5. 07 

CONTINUE 

URITE(21 , 10007 

DO 500 1=1 ,5 

WRITE(21,998) PAR(I7 

CONTINUE 

DO 600 1=1,5 

KCCI7=PAR(I) 

CONTINUE 

CADE RUNGACKC,YC,XC7 
RSS=^FUNK(YC,YM7 
WRITEC21 , 1100) 

WRtTE(21,*7RSS 
WRITR(21 ,997)ITER 

FORMAT C IX, ’THE MATRIX FROM SUBROUTINE AMOEBA IS’ 7 
FORMATCIX, ’THE NUMBER OF ITERATION IS=',I37 
FORMATC1X,1PR10,3) 

FORMATCIX, ’THE VALUE OF PARAMETERS IS') 

FORMATCIX, 'THE RESIDUAL SUM OF SQUARES IS’7 

FORMATCIX, 'THE VALUE FROM RUNGA CALLED IN MAIM PROGRAM’) 

STOP 

END 

**#****-********* + *******A*********A*****A*****’********* 

SUBPROGRAM TO CALCULATE FUNCTION 
FUNCTION FUNKCYC,YM) 

DIMENSION YC(20,3),Yn(20,3) 

SUM-^O 

DO 20 I«l,20 
DO 30 J=l,3 

SUM^SUM+CYMCI , J)-'YCCI , J77**2 

CONTINUE 

CONTINUE 

FUNK^SUM 

RETURN 

END 

**** 1 ^*******^********** ****** **■*'********'** 1 ** ****** 
SUBPROGRAM TO SOLVE SIMULTANEOUS EQUATIONS 
SUBROUTINE RUNGACXC , YC , XC) 

DIMENSION yiC3) ,FC2),YOC2) ,T(4,2) ,yCC20,3),XC(20) 

REAL KCC5) 
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PBO-0. 004 
11== 1 

00 5 S-1.4 
P00=0. 208 
pno-=crROHO 002") 
n=o.o 

XI~-0 

YUt'l-O 

yj c2^-o 

YlO’j^O 
DO 10 I-M,IU-4 
300 DO 100 W=-l,2 

XO=-XJ 

YO(l')-Yia) 
yO(2)«Yl (2) 

CALI, DERlV(KC,YI,F,P00,rB0) 

DO 20 J = 

Tfl, 

20 CONTINUE 

YI(1)^Y0(1)-KH*T(1 , l)/2. O') 

YIC2)'-Y0(2) + (II*TC1 ,2)/2. 0) 

CALL DERIVCKC,Yt ,F,P00,PB0') 

DO s30 J-1 , 2 
Tr2, J)=F(J) 

30 CONTINUE 

XJ‘=X0^ (H/2.0^ 

Yi n :)^yoci)-KH*'i'C2 , i)/2 . 0) 

YI(2)=YO(2) + CH*T(2 , ?.)/2. 0) 

CALL nERIVCKr,YI ,F,P00,PB0) 

DO 40 J~l,2' 

T(3,J)-F(J) 

40 comimv. 

XI«X04H 

YI (3 )•-=YOC■l)^ (H*T(3, 1)) 

YIC2)-=Y0C2) + ('H*T('3, 2)) 

CALL DERIVCKC,YI,F,POO,rBO') 

DO 23 

TC4, J)-F(J) 

23 CONTINUE 
DO 24 

YI(a)=YOCJ)H (T(l, J) + 2*T(2,J)4 2*TC3, J)+T(4,J'))*CH/6.0) 

24 CONTINUE 
YU3)=-YH1)'YI(2) 

100 CONTINUE 

IFfXI.EQ.l.^^GO TO 300 
DO 25 J'-=l,3 
YCCI,J)*=YI(a) 

25 CONTINUE 
XCa)’-=XI 

10 CONTINUE 

n«in>5 

5 CONTINUE 

RETURN 
END 


n o 
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C 


1 


11 


9 


2 


'f(it***icicic‘k‘kifki(-kicic'k-i:icic-k'ki(kkick-klcicic1c-kici(-kif'k1iicic*i(iii(ic'kict-kt'k-kic-k-k 

SUBPROGRAM TO CALCULATE DERIV 
SUBROUTINE DERIVE KC , YI , F , POO , PBO) 

REAL KC(5),F(2),YI(3) 

PB^PBO*(1.0-YI(1)) 

PMA=PBO*YI(2) 

PO-'PCO-a .5*PBO*YIC2)-5.a2*PRO* (yi(l)-YIC2)) 

F(l) = (KC(l)*POi-KC(3)*PB)/(l + (KC(l)*PO/KC(2)*PB) »- 
1 CKC(3)*PB/KCC5)*PO) + (KCC4)*P11A/KC(5)*PO) ) 

FC2) = CKC(l)*PO-KC(a)*PnA)/(l-hCKCC2 :f*PO/KC(2)*FB}-i 
1 fKC(3)*PB/KCC5)*PO)+(KCC4)*PnA/KC(5)*PO)) 

RETURN 

END 

*********************************************************** 
SUBROUTINE AMOEBACP , Y ,MP , NP ,NDIM, FTOL , ITER,YM:) 

PARAMETER CNMAX = 20 , ALPHA=1 . 0 , BETA^O . 5 ,GAMMA^2 . 0 , UMAX =^500 ) 
REAL P(MP ,NP) , YCMP) , PRCNMAX) , PRRCNMAX ) 

1 ,PBAR(NMAX) , YC(20,3) , YM( 20 , 3 ) , XC(20 ) , S(5 ) 

MPTS=NDIM+1 

ITER=0 

ILO^^l 

IF(Y(1) .GT,Y(2 ))THEN 
IHI = 1 
IWHI=-2 
ELSE 
IHI=^2 
INHI=] 

ENDIF 

DO 11 I=1,MPTS 

IFCY(I) .LT.Y(ILO)) ILO-^I 
IFCYCr).GT.ycIHI)lTHEN 
INHI^IHI 
IHI = I 

ELSE IFCYCD .GT.Y(INHI) )THEN 
IFCI.NE.IHI) INHI^I 
ENDIF 
CONTINUE 

RTOL=^2.*ABS(Y(rHI)-'YCILO))/(ABS(y(IHI))^ABS(y(ILO))) 

IF(RTOL.LT.FTOL)RETURN 

IFCITER.EO. ITMAX) RETURN 

ITER=ITER+1 

IF(ITER.EQ.l)GO TO 10 
DO 8 1=1,6 
DO 9 J=l,5 
S(J)"-P(I,J) 

CONTINUE 

CALL RUNGA (S,YC,XC) 

DO 90 K-1,20 
IF(YC(K,1) .GT.1.0)THEN 
URITE(21,2) 

FORMATCIX, 'CONVERSION EXCEEDS UNITY’') 

RETURN 

ELSE 

ENDIF 
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90 CONTINUE 

8 CONTINUE 

iO DO 12 a=l ,NDIM 
PBAR( J)=0. 

12 CONTINUE 

DO 14 1=1,HPTS 
IFCI.NE.IHDTHEN 
DO 13 J=1,NDIN 

PBARCJ)=*PBARCJ) + PCI , J) 

13 CONTINUE 
ENDIF 

14 CONTINUE 

DO 15 J==1,NDIM 

PBAR C J ) =PBAR ( J ) /ND IK 

PR(J)=(I .+ALPHA)*PBAR(J)-ALPHA*PCIHI , J) 

15 CONTINUE 

C YPR-^FUNK(PR) 

CALL RUNGA(PR,YC,XC) 

YPR^FUNKCYC, YM) 

IF(YPR.LE.Y(ILO))THEN 
DO 16 J=l,NDin 

PRRCJ)^GAMMA*PRf J) + (l ,‘-GAMMA)*PBAR( J) 

16 CONTINUE 

C YPRR-=FUNK(PRR) 

CALL RUNGA(PRR,YC,XC) 

YPRR=FUNK(YC,YM) 

IF(YPRR.LT.Y(ILO))THEN 
DO 17 J^I,NDIM 
P(IHI, J)=PRR(J) 

17 CONTINUE 
YaHI)-YPRR 

ELSE 

DO 18 a=l,NDIM 
PCIHI, J)=PR(J) 

18 CONTINUE 
YCIHD^^yPR 

ENDIF 

ELSE IFCYPR.GE,Y(INHI))THBN 
IF(YPR.LT.YCIHI))THEN 
DO 19 J^'1,NDIM 
P(IHI, J)*-PR(J) 

19 CONTINUE 
YCIHI)«YPR 

ENDIF 

DO 21 J==1,NDIM 

PRR(J)-BETA*P(IHI , J) + (:3 . -BETA) * PBAR( J ) 

21 CONTINUE 

C YPRR’-FUNK(PRR) 

CALL RUNGACPRR,YC,XC) 

YPRR^-FUNKCYC.YM) 

IF(yPRR.LT.y(IHI))THEN 
DO 22 J=-1,NDIM 
P(IHI,J)=^PRR(J) 

22 CONTINUE 
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YCIHI)=YPRR 

ELSE 

DO 24 1=1 ,MPTS 
IF(I.NE.ILO)THEN 
DO 23 J = 1,NDII1 

PR(J)^0.5*(PCI , J)+P(IL0, J)) 

PCI , J)'-=PR(J) 

23 COMTIWOE 

C YfI)=FUNK(PR) 

CALL RUNGACPR, YC,XC) 

YCI)=FUNK(YC,YM) 

END IF 

24 CONTINUE 
ENDIF 

ELSE 

DO 25 J=l,NDin 
P(IHI,J)=PR(J) 

25 CONTINUE 
y(IHO=YPR 

ENDIF 
GO TO 1 

100 FORMAT (IX, 'THE ITERATION NUMBER^ 13) 

200 FORMATCIX, 'THE VALUES FROM RUNGA CORRESPONDING TO PRl ’ ) 

300 FORMATCIX, 'THE VALUES FROM RUNGA CORRESPONDING TO PRRl * ) 

400 FORMATCIX, 'THE VALUES FROIi RUNGA CORRESPONDING TO PR2’) 

500 FORMATCIX, 'THE VALUES FROM RUNGA CORRESPONDING TO PRR2*) 

END 

*************5lic*********-**A*****************************<r**** 
r A*********************************************************** 
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